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Abstract:  Efficient  modeling  and  computation 
of  the  nonlinear  interaction  of  fluid  with  a  solid 
undergoing  nonlinear  deformation  has  remained 
a  challenging  problem  in  computational  science 
and  engineering.  Direct  numerical  simulation  of 
the  non-linear  equations,  governing  even  the 
most  simplified  fluid-structure  interaction  model 
depends  on  the  convergence  of  iterative  solvers 
which  in  turn  relies  heavily  on  the  properties  of 
the  coupled  system.  The  purpose  of  this  work  is 
to  model  and  simulate  multi-physics  applications 
that  involve  fluid-structure  interaction  using  a 
distributed  multilevel  algorithm  with  finite 
elements.  The  proposed  algorithm  is  tested  using 
COMSOL  which  offers  the  flexibility  and 
efficiency  to  study  coupled  problems  involving 
fluid-structure  interaction.  Numerical  results  for 
some  benchmark  fluid-structure  interactions  are 
presented  that  validate  the  proposed 
computational  methodology  for  solving  coupled 
problems  involving  fluid-structure  interaction  is 
reliable  and  robust. 

Keywords:  Fluid-Structure  Interaction;  Multi¬ 
physics;  Micro-Air  Vehicle;  Coupled  Systems. 

1.  Introduction 

The  efficient  solution  methodology  to 
complex  multi-physics  problems  involving  fluid- 
structure  interactions  (FSI)  is  a  challenging 
problem  in  computational  sciences.  Such 
problems  require  studying  complex  nonlinear 
interactions  between  the  structural  deformation 
and  the  flow-field  that  often  arise  in  applications 
such  as  blood-flow  interaction  with  an  arterial 
wall  or  computational  aero-elasticity  of  flexible 
micro-air  vehicles. 

In  the  last  two  decades,  domain 
decomposition  techniques  have  become 
increasingly  popular  for  obtaining  fast  and 
accurate  solutions  of  problems  involving  coupled 
processes  [1,  2,  3].  These  viable  domain 
decomposition  techniques  have  been  shown  to  be 
stable  mathematically  and  have  been 
successfully  applied  to  a  variety  of  engineering 


applications  [4,  5,  6].  The  basic  idea  is  to  replace 
the  strong  continuity  condition  at  the  interfaces 
between  the  different  sub-domains  by  a  weaker 
one  to  solve  the  problem  in  a  coupled  fashion. 
The  purpose  of  this  paper  is  to  develop  a  coupled 
FSI  algorithm  and  implement  the  algorithm 
using  COMSOL  to  some  benchmark  FSI 
problems. 

In  section  2,  we  present  the  formulation  of  a 
one -dimensional  FSI  problem  using  an  arbitrary 
Lagrangian  Eulerian  formulation  and  employ  the 
finite  element  method  via  COMSOL  for  solving 
the  coupled  problem.  In  section  3,  we  present  an 
optimal  control  formulation  of  the  FSI  problem 
for  the  1-D  problem  presented  in  section  2.  In 
section  4,  we  consider  a  three  dimensional  FSI 
problem  with  application  to  micro-air  vehicles. 


2.  A  Coupled  One-Dimensional  FSI 
Model  Problem 


For  simplicity  of  presentation,  we  first 
develop  the  model  for  a  one-dimensional  FSI 
problem  that  involves  a  structural  domain 
interacting  with  a  fluid  medium.  The  model  is  set 
up  so  that  initially  the  fluid  domain  occupies  the 
interval  (-1,  0)  and  the  elastic  structure  occupies 
the  interval  (0,  1).  As  the  fluid  flow  deforms  the 
adjacent  solid,  we  allow  the  movement  of  the 
interface  to  depend  on  the  velocity  of  the  fluid. 
This  is  illustrated  in  figure  1  below. 
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Figure  1:  Undeformed  and  deformed 
computational  domains 

We  denote  by  y(t)  the  position  of  this 
interface  at  any  positive  time  t.  For  all  values  of 
x  in  the  interval  (-1,  y(t ) ),  we  model  the  fluid 
velocity  v  and  the  pressure  p  using  a 
generalization  of  models  employed  by  [7]  using: 

v,-avB  +  (l  +  p)vvx  +  £Px=ff 
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In  this  model,  OC  =  C((v)  >  0  is  a  parameter 
depending  on  the  kinematic  viscosity  V  of  the 
fluid.  The  constants  /?  E  [0,0.5]  and 
£  >  0  vary  depending  on  the  material  properties 
of  the  fluid.  Here  j ;  is  the  external  force  on  the 

fluid.  In  this  work  we  consider  the  specific  case 
of  the  Burger’s  equation  using  the  parameter 
choices,  a  =  V,  (3  =  0.5,  £  =  0  .  It  is  well 
known  that  this  choice  of  values  yields  stability 
of  the  numerical  methods  employed  for  the 
problem  whenever  no  pressure  term  is  included. 

The  fluid  model  is  coupled  with  an  elastic 
model  that  represents  the  structure.  In  particular, 
we  consider  the  wave  equation  that  models  the 
solid  displacement  u  of  any  point  in  the  adjacent 
structure  from  its  initial  position  given  by: 

Utt~/iUxx  =fs 

Here  ju  >  0  and  fs  is  the  external  force  on 

the  structure.  Also,  the  position  of  the  interface 
between  the  two  sub-domains  must  satisfy  the 
movement  at  the  interface  for  all  times. 

At  the  interface  between  the  fluid  and 
structure,  we  enforce  continuity  of  the  fluid 
velocities  and  the  action-reaction  principle: 


v(-«vxc+(1-5v-wK=//' 

For  the  boundary  conditions,  we  let  the  fluid 
velocity  v(-l,  t)=sin(0. 1  71 1)  and  the  structural 
displacement  u(  1,  t)  =  0  for  all  times.  For  initial 
conditions,  both  the  fluid  and  the  structural 
domain  are  at  rest  initially. 

The  FS1  coupled  system  described  was 
implemented  in  COMSOL  and  the  results  of  the 
finite  element  implementation  of  the  model  is 
summarized  in  figure  2  that  displays  the  plot  of 
the  fluid  velocity  v  and  the  structural  velocity  u, 
over  the  time  period  from  0  to  1 . 


SOW  dupliccment  velocity 


Figure  2:  FSI  simulation  using  COMSOL 
for  time  t=0  to  1 . 


3.  Coupled  FSI  problem  with  Control 


a  0(0,0  =  /*ux(Q,t) 

In  order  to  account  for  the  changing  nature  of 
the  fluid,  we  consider  the  arbitrary  Lagrangian- 
Eulerian  (ALE)  formulation  [8].  This  will  allow 
for  a  dynamic  computational  gird  that  avoids 
extreme  mesh  distortion  near  the  interface.  To  do 
this  one  can  move  the  numerical  grid 
independently  of  the  fluid  velocity  on  the  fluid 
domain.  Defining  the  grid  velocity  as: 

X 

w  =  — 


Note  that t)  =  y{t),  w(0,f)  =  0.  Thus 
the  grid  velocity  is  consistent  with  the  velocity  of 
the  fluid  at  the  endpoints  of  the  fluid  domain. 
Additionally,  we  assume  E  [— 1,1]  for  all 
time.  The  ALE  form  of  the  fluid  equation  that  we 
then  solve  is: 


A  related  aspect  that  we  consider  in  this  work 
is  investigating  distributed  control  for  FSI 
problems.  In  particular,  one  can  study  the 
coupled  FSI  problem  using  an  optimal  control 
formulation  to  predict  the  distributed  control  that 
corresponds  to  a  prescribed  velocity  Vand 
displacement  data  U  that  satisfies  the  boundary 
conditions.  Specifically,  we  want  develop  a 
model  that  predicts  the  force  that  results  in  the 
minimization  of  the  error  in  the  fluid  velocity 
|  V  —  V  |  and  the  solid  displacement  |  U  —  U  \ . 
For  simplicity,  let  us  consider  the  one¬ 
dimensional  model  problem  presented  in  section 
2  and  extend  it  to  include  control  aspects. 

Towards  this  end,  let  us  consider  the  related 
cost  functional  for  the  associated  non-linear  FSI 
problem  (note  that  one  can  similarly  consider  the 
linear  FSI  problem  by  dropping  the 
corresponding  linear  terms  which  is  not 
described  here  )  given  by: 


M  =  \\[ \b-vf+\afff 

o-iVz  z  y 

+  jjfy(w-«)2+y«s/; 

O-lV^  ^  J 

T  0 

+  J  j  (KpfV,  ~  Mfvxx  +  1  -5vvx  -  ffj) 
0  -1 
T  0 

+  J  j  (g(PsU«  ~  MsUxx  ~  f,)) 


0  -1 


where  l  and  g  are  the  Lagrange  multipliers 
corresponding  to  the  fluid  velocity  and  the  solid 
displacement  respectively.  Moreover,  we  also 
impose  the  continuity  of  the  velocities  over  the 
coupled  domain. 

Proceeding  using  the  standard  optimal 
control  approach  of  minimizing  the  cost  function 
by  taking  the  variations  yields  the  following 
auxiliary  system  of  governing  equations: 

In  the  fluid  domain  —  1  <  X  <  0  : 


P/Vt-^fv„  +  l-5wx - =  0 

af 

~  Pfvt  -  Mfl„  -  !-5  vlx  +  v  -  v  =  0 


In  the  solid  domain  0  <  X  <  1 : 

Psvt  -  Psuxx  -  —  =  0 
PsSa  -  PsSxx  +u-u  =  0 


For  numerical  experiments,  we  consider  a 
simple  solid  displacement  profile  and  the  fluid 
velocity  profile  given  by: 

u  =0.5.r(x2  -  l)f2 
v  =x(x2  -  l)t 

The  problem  was  implemented  in  COMSOL  to 
yield  the  following  velocity  and  displacement 
profiles  in  the  respective  domains.  This  is 
illustrated  in  Figures  3  and  4  where  the 
prescribed  solution  is  plotted  against  both  the 
solution  to  the  both  the  linear  and  non-linear 
control  problem. 
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Figure  3:  Displacement  profile  comparing  the 
prescribed  solid  displacement  with  linear  and 
non-linear  control  models. 


Figure  4:  Fluid-velocity  profile  comparing  the 
prescribed  velocity  with  linear  and 
non-linear  control  models 

4.  A  Coupled  multi-dimensional  FSI  problem 
with  applications  to  MAVs 

A  Micro  Air  Vehicle  (MAV)  is  a  type  of  radio- 
controlled  miniature  aircraft  that  can  fly  at  very 
low  speeds.  Due  to  the  complexities  of  the  wing 
structure  of  a  MAV,  a  computational  model  of 
the  aircraft  wing  requires  a  combination  of  many 
structural  elements  interacting  with  external 
fluid. 

The  wing  typically  consists  of  a  flexible 
membrane  material  braced  with  a  leading  edge 
spar  and  chordwise  battens  (see  Figure  5).  The 
structural  model  must  combine  the  model  of  the 
membrane  material  together  with  the  model  of 
the  rigid  battens.  Most  current  models  treat  the 
battens  as  large-density  membrane  elements. 
Modeling  this  coupled  two-dimensional 
structural  model  interacting  with  a  three- 
dimensional  fluid  makes  the  problem  very 
challenging. 


Figure  5:  A  simple  model  of  a  flexible  MAV  wing 

In  this  work,  we  attempt  to  come  up  with  a 
mathematical  model  that  can  provide  insight  into 
the  dynamics  of  MAVs.  The  model  presented 
herein  is  simplistic;  however,  one  may  extend 
this  to  accommodate  other  features. 

In  order  to  get  an  insight  into  the  modeling 
and  dynamics  of  a  MAV,  let  us  consider  a 
cylindrical  computational  domain.  In  this  domain 
we  will  assume  that  the  fluid  satisfies  a  potential 
equation  and  that  a  two-dimensional  structural 
model  (that  will  represent  the  MAV)  is  a  part  of 
one  of  the  circular  surfaces.  This  latter  surface 
will  represent  the  outflow  boundary  of  the 
computational  domain  which  is  illustrated 
below: 


Figure  6:  Computational  domain  for  MAVs 


Let  the  computational  domain  be  partitioned  into 
three  sub-domains  Q(. ,  (=1,2,3.  Let  Q , 

represent  the  cylinder  in  the  computational 
domain  where  the  following  governing  equations 
hold: 


A$)  =  0  in  L23 

V(f>-n- o  on  r;' 

Here  corresponds  to  the  lateral  surface  of 

the  cylinder  where  Neumann  boundary 
conditions  are  prescribed. 

The  outflow  part  of  the  computational 
domain  consist  of  the  following  sub-domains: 

T°  corresponds  to  the  outflow  region  that  is  not 

a  part  of  the  structural  domain;  Q2  corresponds 
to  the  structural  domain  that  involves  the  three 
battens;  O,  corresponds  to  the  structure  (shaded 
grey)  that  does  not  involve  the  three  battens.  We 
will  assume  the  following  adsorption  condition: 

V (j)-n=-a(j)t  on  T° 

where  a  is  a  constant.  Also,  IA  corresponds  to 
the  inflow  surface  where  we  prescribe: 

V^-/7=-0.1  +  0.025sin(2^f)  on  T‘f 

For  the  structural  model  of  the  MAV  that  is 
modeled  via  the  sub-domains  Qj  and  Q,  we 
consider  the  following  governing  membrane 
equations  for  the  deflection  of  the  membrane  w : 

Po  wtt  ~  E0  Aw  =  /  in  Q, 

(  a,  +  pl  )wll-Ea  Aw  +  Et  v„,  =  -  pf  fa  in  Q2 

v=w  +  s Av  in  Q,uQ 

where  E0 ,  E{  are  the  constants  corresponding  to 
the  elastic  modulus  and  second  moment  of  area 
for  each  of  the  sub-domains  Q[  and  Q2 
respectively.  Also,  PG,P\  are  the  respective 

densities  of  the  membrane  and  the  battens.  It 
must  also  be  pointed  out  that  half  of  the  MAV 
edge  was  kept  rigid  to  reflect  the  leading  edge 
spar.  The  two  systems  are  also  coupled  through 
the  continuity  of  the  velocities: 

V(f)-n  =  wt  on  Q, 

The  fully  coupled  system  described  herein  was 
modeled  and  solved  in  COMSOL  and  the  results 
for  the  membrane  deflection  are  shown  in 
Figures  7  and  8. 


Time=3  Max:  2.087 


Min:  -4.466e-4 

Figure  7:  Three-dimensional  MAV 
membrane  wing  deflection 

Time=3  Max:  2.087 


Figure  8:  Two-dimensional  MAV 
membrane  wing  deflection 

5.  Conclusions 

In  this  work,  a  coupled  computational 
methodology  to  solve  problems  that  involve 
fluid-structure  interaction  has  been  presented  for 
various  benchmark  problems.  The  problems 
considered  in  this  work  included  a  one 
dimension  problem  coupling  fluid  and  structure 
with  and  without  control  and  an  application 
problem  in  three  dimensions  involving  MAVs. 
The  one  dimensional  problems  provide  a  great 
insight  into  the  nature  of  the  coupled  behavior  of 
the  interaction  between  the  fluid  velocity  and  the 
structural  displacement.  The  importance  of  the 
non-linear  term  in  the  fluid  equations  was 
illustrated  in  the  control  problem  that  helped 
decrease  the  error  between  the  prescribed  and 
computed  solutions. 
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Flow-structure  interaction  (FSI) 
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Fluid-Structure  Interaction  (ID) 
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Coupled  FSI  with  Control 
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Auxiliary  system  of  PDEs 


In  the  fluid  domain: 


i 


PfUt~PfUxx  +  L5  UUx - =  0 


a 


f 


pfut  —  juflxx—\  .5  ulx  +  u  —  u  —  0 


In  the  solid  domain: 


8 


PsU,  -  - =  0 


a 


Ps8tt-M.s8xx+d-d  =  0 


Fluid  Velocity  Profile 


Fluid-velocity  Profile 

d  =  0.5x(x2  -  l)^2 
u  =  x{x2  -  \)t 


Displacement  Profile 

d  =  0.5x(x2  -  l)^2 
u  =  x{x2  -  \)t 


MAV:  Membrane  Wing  Deflection 


Computational  modeling  of  highly  flexible  membrane  wings  for  MAVs 

Ferguson,  L.,  Aulisa,  E.,  Seshaiyer  P.  and  Gordnier  R.,  (AIAA  2006-1 661) 

Computational  modeling  of  coupled  membrane-beam  flexible  wings  for  MAVs 

L.  Ferguson,  P.  Seshaiyer,  R.  Gordnier,  P.  Attar  (AIAA  2007-1787) 

Nonlinear  Models  for  Biologically-  Inspired  Elastic  Membrane  Wings 

E.  Swim  and  P.  Seshaiyer  (AIAA-2008-2008) 


Experimental  Challenges  in 

MAV  Design 

Small  size 

High  surface-to-volume  ratio 

Constrained  weight  and  volume  limitations 

Low  Reynolds  number  regime 

Low  aspect  ratio  fixed  to  rotary  to  flapping 
wings 

Longer  flight  time 

Better  range-payload  performance 


Computational  Challenges  in 

MAV  Design 


Structural 
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Dynamics 
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Model  of  a  flexible  MAV  wing 


Membrane 

Material 


Leading  Edge 
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Computational  Model  for  MAV 


A^=0  in  Q3 
V^-n  = 0  on  r/ 


A  M;„  ~  E0  Aw  =  f  in  Q, 

(A+A)M;a-£oAM;  +  £,ivXy=-/7/^  in  Q2 

v  =  wyy  +  oAv  zn  QiUQ2 


V cj)-n—  —  a (j){  on  T° 

V(f)-n=—  0.1  +  0.025 sin(2^)  on 
V(f)-n  =  wt  on  Qj 


Membrane  Wing  Deflection 

Time=3  Max:  2.087 


Min:  -4.466e-4 


Membrane  Wing  Deflection 


Time=3 
Surface:  w 
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